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Abstract 

The cumulants and moments of the log of the non-central chi-square 
distribution are derived. For example, the expected log of a chi-square 
random variable with v degrees of freedom is log 2 + 1 /; (|). Applications 
to modeling probability distributions are discussed. 


1 Introduction 

The Edgeworth and Cornish-Fisher expansions allow one to approximate the 
density, distribution, and quantile functions of probability distributions whose 
cumulants are known. It is often remarked that these expansions are inaccurate 
for one-sided and highly skewed probability distributions. [2, 3] 

For example, consider a random variable which is the product of several 
chi-square random variables. The raw moment of a chi-square random 
variable with v degrees of freedom is 2^r {k v/2) /T (v/2). The raw moments 
of the product of multiple independent chi-squares is simply the product of 
the moments. The moments can then be converted to the cumulants. The 
cumulants are then used in the Edgeworth expansion to approximate the density. 
For even a small number of factors, the Edgeworth expansion is inaccurate, 
sometimes yielding negative estimates, as illustrated in Figure 1. 

# moments of chi-square 

chisq.moments <- function(df, order.max = 3) { 
orders <- 1 :order.max 

mu <- exp (orders * log (2) + Igamma (orders + 

(df/2)) - lgamma(df /2) ) 
return (mu) 

} 

# moments of product of chi-squares 

prodchisq.moments <- function(df s , order.max = 3) { 
mu <- Reduce ("*", sapply(dfs, chisq.moments, 
order.max = order.max, simplify = FALSE)) 
return (mu) 

} 

rprodchisq <- function(n, dfs) { 
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X <- Reducesupply (dfs, function(nu) { 
rchisqCn, df = nu) 

}, simplify = FALSE)) 
return (X) 

} 

require (PDQutils) 

dprodchisq <- function(x, dfs, log = FALSE) { 

kappa <- PDQutils: :moment2cumulant(prodchisq.moments (dfs, 
order.max = 4)) 

pdf <- dapx_edgeworth(x, kappa, support = c(0, 

Inf), log = log) 
return (pdf) 

} 

require (ggplot2) 

test_dens <- function(dpqr , nobs, ...) { 
rv <- sort (dpqr$r(nobs, ...)) 
data <- data.frame (draws = rv) 
text.size <- 6 # sigh 

# I widthOpthttp ://stackoverflow. com/a/5688125/16J,.611 
pi <- qplot(rv, geom = "blank") + 
geom_line(aes (y = ..density.., 

colour = "Empirical"), stat = "density") + 
stat_function (fun = function(x) { 
dpqr$d(x, ...) 

}, aes (colour = "Theoretical")) + 
geom_histogrcLm(aes (y = ..density..), 

alpha = 0.3) + scale_colour_manual (name = "Density", 
values = cC'red", "blue")) + 

theme (text = element_text (size = text.size)) + 
labs (title = "Density (tests dfunc)") 
print (pi) 

} 

dfs <- c(40, 30, 50, 20, 10) 

test_dens(list (r = rprodchisq, d = dprodchisq), 
nobs = 2^14, dfs) 


2 Moments of the log chi-square distribution 

One possible fix to this problem is to use the Edgeworth expansions to approx¬ 
imate the density of the log of the chi-square, or the weighted sum of logs of 
chi-square variates. This latter approach would allow one to model the dou¬ 
bly non-central F distribution, for example, as well as products of arbitrary 
chi-squares to different powers. 

Let X ^ (t^) be a chi-square variate with v degrees of freedom. Let y = 

logx. Consider the cumulant generating function of It is 

K{t) =logE [e^^] =logE [x^]. 
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Density (tests dfunc) 



Figure 1: The Edgeworth expansion is inaccurate for the product of chi-squares 
variates. 


The moments of the central chi-square are known [6], yielding the expression 
K {t) = Hog2 + logr (^ + *) - logT 

The cumulant of y is the derivative evaluated at t = 0. The derivative of 
the log of the Gamma function is the ‘psi’ function, (•): while its derivatives 
are the ‘polygamma’ functions, (•). [1, 6.3, 6.4] Letting kj be the raw 
cumulant of y, we have 


log2 + '^(|), 

(I) , 


if j = 1, 

if j > 1. 


( 1 ) 


The moments can then be computed from the cumulants via the usual for¬ 
mula: 



Note that the first raw cumulant equals the first moment. That is. 


E[y] = log2 + ?/) (^1) . 


(3) 


3 Moments of the log non-central chi-square dis¬ 
tribution 

To compute the moments of the log of the non-central chi-square, the crucial 
observation is that the density of the non-central chi-square can be expressed as 
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a Poisson mixture of central chi-squares. [6] That is, if f ^2 (x;v) is the density of 
the central chi-square distribution with v degrees of freedom, and /^ 2 / (x; A, v) 
is the density of the non-central chi-square with v degrees of freedom and non¬ 
centrality parameter A, then 

oo /Ay’ 

f^ 2 , {x; X,v)=Y^ {x; v + 2j). (4) 

i=o 

Using the change of variables formula, if y is the log of a non-central chi- 
square variate with v degrees of freedom and non-centrality parameter A, then 
the density of y follows a similar relationship: 

oo /Ay’ 

/y (y; A, v) = eyf^ 2 , (e^; A, v)=e'^Y^ (e^; v + 2j). (5) 

j=o 


Because the uncentered moments are defined in terms of an integral, which 
is a linear operator, the moments of y can be expressed as a similar sum: 


E 


oo /Ay nOO 

= E / e^/x 2 (e^; v + 2j) z'^ dz, 

n—n J—oo 


j=0 
oo 

= E^ 

i=o 


-A/2 


(6) 


J 


I k,v-\-2j 5 


where y'^ 2 j moment of the log of the central chi-square distribution 

with V 2j degrees of freedom. 

For the particular case of /c = 1, because the first cumulant is the first 
moment, via Equation 1, we have 

00 /Ay’ 

E[y] = log2 + ^e-V2l2i_^(j + ^/2), (7) 

j=o 

Note that this does not seem to match the equations given by Moser, even in 
the simple case A = 0. [4] (It appears that Moser’s result is missing a summand 
of log 2.) It is easy to check this formula via Monte Carlo simulations, as below. 
The results of these experiments, reported in Table 1 indicate that the equations 
are indeed accurate. 


require (PDQutils) 

# cumulants of the log of the central 

# chi-square 

lc_cumuls <- function(df, order.max = 3, 
orders = c (1 :order.max)) { 
kappa <- psigamma(df /2 , deriv = orders - 
1 ) 

kappa[1] <- kappa[1] + log (2) 
return (kappa) 

} 

# moments of the log of the central 
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# chi-square 

Ic.moments <- function (df , order.max = 3, 
orders = c (1 :order.max)) { 
kappa <- lc_cumuls (df, orders = orders) 
mu <- PDQutils:: cumulant2moment (kappa) 
return (mu) 

} 

# moments of the log of the 

# non-central chi-square 
lnc_moments <- function (df, ncp = 0, 

order.max = 3, orders = c (1 :order.max)) { 
stopifnot (ncp >= 0) 
if (ncp > 0) { 

hancp <- ncp/2 

# should he smarter about 0:100 here. 
allmu <- sapply (0 : 100 , function (iv) { 
exp (-hancp + iv * log (hancp) - 

Ifactorial (iv)) * lc_moments (df + 
2 * iv, orders = orders) 

}, simplify = FALSE) 
mu <- Reduce(" + '', allmu) 

} else { 

mu <- lc_moments (df = df, orders = orders) 

} 

return (mu) 

} 

set. seed(1231591) 
df <- 50 
ncp <- 1.5 
nsim <- le+06 

X <- rchisq(nsim, df = df, ncp = ncp) 
y <- log(x) 

nord <- 6 

empirical.mu <- sapply (1 :nord, function(k) { 
mean(y"k) 

}) 

theoretical.mu <- lnc_moments (df = df, 
ncp = ncp, order.max = nord) 


4 Using the Moments 

While the moments computation is perhaps of theoretical interest, the nominal 
impetus for this work was more accurate simulation of the density of prod¬ 
ucts of non-central chi-squares taken to powers. Here we first approximate the 
density of the log of such a distribution, using additivity of cumulants, via an 
Edgeworth expansion, then use change of variables to recover the density of the 
product of chi-squares. The resultant density estimator, shown in Figure 2, is 
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order 

empirical 

theoretical 

1 

3.92 

3.92 

2 

15.42 

15.42 

3 

60.79 

60.78 

4 

240.26 

240.22 

5 

951.97 

951.78 

6 

3781.22 

3780.36 


Table 1: Empirical and theoretical moments up to order 6 are shown for 10^ 
draws from the log of a non-central chi-square distribution with 50 degrees of 
freedom and non-centrality parameter 1.5. 


an improvement, at least under an ‘eyeball test.’ 

# cumulants of the log of the 

# non-central chi-square 

liic_cumuls <- function(df , ncp = 0, order.max = 3, 
orders = c (1 :order.max)) { 

mu <- liic_moments (df, ncp, orders = orders) 
kappa <- PDQutils:: moment2cumulant (mu) 
return (kappa) 

} 

# compute the cumulants of the 

# sumlogchisq distribution. 
sumlogchisq_cumuls <- function(wts , df, 

ncp = 0, order.max = 3) { 
subkappa <- mapply (function (w, dd, 
nn) { 

(w" (1 :order.max)) * lnc_cumuls (df = dd, 
ncp = nn, order.max = order.max) 

}, wts, df, ncp, SIMPLIFY = FALSE) 
kappa <- Reduce("+", subkappa) 
return (kappa) 

} 

dsumlogchisq <- function(x, wts, df, 

ncp = 0, log = FALSE, order.max =6) { 
kappa <- sumIogchisq_cumuIs (wts, 

df, ncp, order.max = order.max) 
retval <- PDQutils: :dapx_edgeworth(x, 
kappa, log = log) 
return (retval) 

} 

# use change of variables: 

dprodchisq2 <- function(x, dfs, log = FALSE) { 
dx <- dsumIogchisq(Iog(x) , wts = 1, 
df = dfs, ncp = 0, log = log) 
if (log) { 

dx <- dx - log(x) 

} else { 

dx <- dx/x 

} 
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Density (tests dfunc) 
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Figure 2: Applying the Edgeworth expansion to the sum of log chi-squares, then 
performing a change of variables, results in a more accurate density estimate. 


return (dx) 

} 

dfs <- c(40, 30, 50, 20, 10) 

test_dens(list (r = rprodchisq, d = dprodchisq2), 
nobs = 2^14, dfs) 
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